APIM - Control - Pushing Separate
Pascal Küng, MSc, Patrick Höhener, MSc, James M. Allen, PhD, Robert Tobias, PhD, Urte Scholz, PhD
2024-11-24
knitr::opts_chunk$set(echo = TRUE)
cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)
source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))
report_function_hash <- digest::digest(summarize_brms)## [1] 1116# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
report_hurdle = TRUE
do_priorsense = FALSE
options(
dplyr.print_max = 100,
brms.backend = 'cmdstan',
brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
brms.file_refit = 'on_change',
#brms.file_refit = 'always',
error = function() {
beepr::beep(sound = 5)
if (shutdown) {
system("shutdown /s /t 180")
quit(save = "no", status = 1)
}
}
, es.use_symbols = TRUE
)
####################### Model parameters #######################
iterations = 2000 # 10'000 per chain to achieve 40'000
warmup = 1000
# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'
################################################################
suffix_nopushing = paste0('_sensitivityNoPushing_', as.character(iterations))
suffix_onlypushing = paste0('_sensitivityOnlyPushing_', as.character(iterations))df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df
df_double <- prepare_data(df, recode_pushing = FALSE, use_mi = use_mi)[[1]]Constructing scales reshaping data (4field) centering data within and between
Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.3468 0.0000 5.0000 2323
Modelling
# For indistinguishable Dyads
model_rows_fixed <- c(
'Intercept',
# '-- WITHIN PERSON MAIN EFFECTS --',
'persuasion_self_cw',
'persuasion_partner_cw',
'pressure_self_cw',
'pressure_partner_cw',
'pushing_self_cw',
'pushing_partner_cw',
'day',
'plan_selfPlan',
'plan_partnerPlan',
'weartime_self_cw',
# '-- BETWEEN PERSON MAIN EFFECTS',
'persuasion_self_cb',
'persuasion_partner_cb',
'pressure_self_cb',
'pressure_partner_cb',
'pushing_self_cb',
'pushing_partner_cb',
'weartime_self_cb'
)
model_rows_fixed_ordinal <- c(
model_rows_fixed[1],
'Intercept[1]',
'Intercept[2]',
'Intercept[3]',
'Intercept[4]',
'Intercept[5]',
model_rows_fixed[2:length(model_rows_fixed)]
)
model_rows_random <- c(
# '--------------',
# '-- RANDOM EFFECTS --',
'sd(Intercept)',
'sd(persuasion_self_cw)',
'sd(persuasion_partner_cw)',
'sd(pressure_self_cw)',
'sd(pressure_partner_cw)',
'sd(pushing_self_cw)',
'sd(pushing_partner_cw)',
# '-- CORRELATION STRUCTURE -- ',
'sigma'
)
model_rows_random_ordinal <- c(model_rows_random,'disc')# For indistinguishable Dyads
model_rownames_fixed <- c(
"Intercept",
# "-- WITHIN PERSON MAIN EFFECTS --",
"Daily persuasion experienced",
"Daily persuasion utilized (partner's view)", # OR partner received
"Daily pressure experienced",
"Daily pressure utilized (partner's view)",
"Daily pushing experienced",
"Daily pushing utilized (partner's view)",
"Day",
"Own Actionplan",
'Partner Actionplan',
"Daily weartime",
# "-- BETWEEN PERSON MAIN EFFECTS",
"Mean persuasion experienced",
"Mean persuasion utilized (partner's view)",
"Mean pressure experienced",
"Mean pressure utilized (partner's view)",
"Mean pushing experienced",
"Mean pushing utilized (partner's view)",
"Mean weartime"
)
model_rownames_fixed_ordinal <- c(
model_rownames_fixed[1],
'Intercept[1]',
'Intercept[2]',
'Intercept[3]',
'Intercept[4]',
'Intercept[5]',
model_rownames_fixed[2:length(model_rownames_fixed)]
)
model_rownames_random <- c(
# '--------------',
# '-- RANDOM EFFECTS --',
'sd(Intercept)',
"sd(Daily persuasion experienced)",
"sd(Daily persuasion utilized (partner's view))", # OR partner received
"sd(Daily pressure experienced)",
"sd(Daily pressure utilized (partner's view))",
"sd(Daily pushing experienced)",
"sd(Daily pushing utilized (partner's view))",
# '-- CORRELATION STRUCTURE -- ',
'sigma'
)
model_rownames_random_ordinal <- c(model_rownames_random,'disc')rows_to_pack <- list(
"Within-Person Effects" = c(2,10),
"Between-Person Effects" = c(11,17),
"Random Effects" = c(18, 24),
"Additional Parameters" = c(25,25)
)
rows_to_pack_ordinal <- list(
"Intercepts" = c(1,6),
"Within-Person Effects" = c(2+5,10+5),
"Between-Person Effects" = c(11+5,17+5),
"Random Effects" = c(18+5, 24+5),
"Additional Parameters" = c(25+5,25+6)
)HURDLE MODELS
# For indistinguishable Dyads
model_rows_fixed_hu <- c(
'Intercept',
'hu_Intercept',
# '-- WITHIN PERSON MAIN EFFECTS --',
'persuasion_self_cw',
'persuasion_partner_cw',
'pressure_self_cw',
'pressure_partner_cw',
'pushing_self_cw',
'pushing_partner_cw',
'day',
'plan_selfPlan',
'plan_partnerPlan',
'weartime_self_cw',
# '-- BETWEEN PERSON MAIN EFFECTS',
'persuasion_self_cb',
'persuasion_partner_cb',
'pressure_self_cb',
'pressure_partner_cb',
'pushing_self_cb',
'pushing_partner_cb',
'weartime_self_cb',
# HURDLE MODEL
# '-- WITHIN PERSON MAIN EFFECTS --',
'hu_persuasion_self_cw',
'hu_persuasion_partner_cw',
'hu_pressure_self_cw',
'hu_pressure_partner_cw',
'hu_pushing_self_cw',
'hu_pushing_partner_cw',
'hu_day',
'hu_plan_selfPlan',
'hu_plan_partnerPlan',
'hu_weartime_self_cw',
# '-- BETWEEN PERSON MAIN EFFECTS',
'hu_persuasion_self_cb',
'hu_persuasion_partner_cb',
'hu_pressure_self_cb',
'hu_pressure_partner_cb',
'hu_pushing_self_cb',
'hu_pushing_partner_cb',
'hu_weartime_self_cb'
)
model_rows_fixed_hu_ordinal <- c(
model_rows_fixed_hu[1:2],
'Intercept[1]',
'Intercept[2]',
'Intercept[3]',
'Intercept[4]',
'Intercept[5]',
model_rows_fixed_hu[3:length(model_rows_fixed_hu)]
)
model_rows_random_hu <- c(
# '--------------',
# '-- RANDOM EFFECTS --',
'sd(Intercept)',
'sd(hu_Intercept)',
'sd(persuasion_self_cw)',
'sd(persuasion_partner_cw)',
'sd(pressure_self_cw)',
'sd(pressure_partner_cw)',
'sd(pushing_self_cw)',
'sd(pushing_partner_cw)',
# HURDLE
'sd(hu_persuasion_self_cw)',
'sd(hu_persuasion_partner_cw)',
'sd(hu_pressure_self_cw)',
'sd(hu_pressure_partner_cw)',
'sd(hu_pushing_self_cw)',
'sd(hu_pushing_partner_cw)',
# '-- CORRELATION STRUCTURE -- ',
'sigma'
)
model_rows_random_hu_ordinal <- c(model_rows_random_hu,'disc')# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
"Intercept",
"Hurdle Intercept",
# "-- WITHIN PERSON MAIN EFFECTS --",
"Daily persuasion experienced",
"Daily persuasion utilized (partner's view)", # OR partner received
"Daily pressure experienced",
"Daily pressure utilized (partner's view)",
"Daily pushing experienced",
"Daily pushing utilized (partner's view)",
"Day",
"Own actionplan",
'Partner actionplan',
"Daily weartime",
# "-- BETWEEN PERSON MAIN EFFECTS",
"Mean persuasion experienced",
"Mean persuasion utilized (partner's view)",
"Mean pressure experienced",
"Mean pressure utilized (partner's view)",
"Mean pushing experienced",
"Mean pushing utilized (partner's view)",
"Mean weartime",
# HURDLE
# "-- WITHIN PERSON MAIN EFFECTS --",
"Hu Daily persuasion experienced",
"Hu Daily persuasion utilized (partner's view)", # OR partner received
"Hu Daily pressure experienced",
"Hu Daily pressure utilized (partner's view)",
"Hu Daily pushing experienced",
"Hu Daily pushing utilized (partner's view)",
"Hu Day",
"Hu Own actionplan",
'Hu Partner actionplan',
"Hu Daily weartime",
# "-- BETWEEN PERSON MAIN EFFECTS",
"Hu Mean persuasion experienced",
"Hu Mean persuasion utilized (partner's view)",
"Hu Mean pressure experienced",
"Hu Mean pressure utilized (partner's view)",
"Hu Mean pushing experienced",
"Hu Mean pushing utilized (partner's view)",
"Hu Mean weartime"
)
model_rownames_fixed_hu_ordinal <- c(
model_rownames_fixed_hu[1:2],
'Intercept[1]',
'Intercept[2]',
'Intercept[3]',
'Intercept[4]',
'Intercept[5]',
model_rownames_fixed_hu[3:length(model_rownames_fixed_hu)]
)
model_rownames_random_hu <- c(
# '--------------',
# '-- RANDOM EFFECTS --',
'sd(Intercept)',
'sd(Hurdle Intercept)',
"sd(Daily persuasion experienced)",
"sd(Daily persuasion utilized (partner's view))", # OR partner received
"sd(Daily pressure experienced)",
"sd(Daily pressure utilized (partner's view))",
"sd(Daily pushing experienced)",
"sd(Daily pushing utilized (partner's view))",
# Hurdle
"sd(Hu Daily persuasion experienced)",
"sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
"sd(Hu Daily pressure experienced)",
"sd(Hu Daily pressure utilized (partner's view))",
"sd(Hu Daily pushing experienced)",
"sd(Hu Daily pushing utilized (partner's view))",
# '-- CORRELATION STRUCTURE -- ',
'sigma'
)
model_rownames_random_hu_ordinal <- c(model_rownames_random_hu,'disc')rows_to_pack_hu <- list(
"Conditional Within-Person Effects" = c(3,12),
"Conditional Between-Person Effects" = c(13,19),
"Hurdle Within-Person Effects" = c(20,29),
"Hurdle Between-Person Effects" = c(30,36),
"Random Effects" = c(37, 50),
"Additional Parameters" = c(51,51)
)
rows_to_pack_hu_ordinal <- list(
"Conditional Within-Person Effects" = c(3+5,12+5),
"Conditional Between-Person Effects" = c(13+5,19+5),
"Hurdle Within-Person Effects" = c(20+5,29+5),
"Hurdle Between-Person Effects" = c(30+5,36+5),
"Random Effects" = c(37+5, 50+5),
"Additional Parameters" = c(51+5,51+6)
)WITHOUT PUSHING
Self-Reported MVPA
Hurdle Lognormal Model
formula <- bf(
pa_sub ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
plan_self + plan_partner +
day +
# Random effects
(1 + persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw | coupleID),
hu = ~ persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
plan_self + plan_partner +
day +
# Random effects
(1 + persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
, brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix_nopushing))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
pa_sub,
stats_to_report = c('pd','ROPE'),
rope_range = c(-0.18, 0.18),
model_rows_fixed = model_rows_fixed_hu,
model_rows_random = model_rows_random_hu,
model_rownames_fixed = model_rownames_fixed_hu,
model_rownames_random = model_rownames_random_hu,
exponentiate = T) ## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub_continuous <- summarize_brms(
pa_sub,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
rope_range = rope_range_continuous,
model_rows_fixed = model_rows_fixed_hu,
model_rows_random = model_rows_random_hu,
model_rownames_fixed = model_rownames_fixed_hu,
model_rownames_random = model_rownames_random_hu,
exponentiate = T) ## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous
columns_to_replace <- c("ROPE", "inside ROPE")
summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <-
summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]
# Print the updated dataframe
summary_pa_sub %>%
print_df(rows_to_pack = rows_to_pack_hu)| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 38.46*** | 2.69 | [33.40, 44.04] | 1.000 | [0.92, 1.08] | 0.000 | 1.003 | 1203 | 2271 |
| Hurdle Intercept | 0.25*** | 0.04 | [ 0.19, 0.35] | 1.000 | [0.84, 1.20] | 0.000 | 1.000 | 1088 | 1791 |
| Conditional Within-Person Effects | |||||||||
| Daily persuasion experienced | 1.02 | 0.03 | [ 0.97, 1.08] | 0.817 | [0.92, 1.08] | 0.973 | 1.002 | 1741 | 2793 |
| Daily persuasion utilized (partner’s view) | 1.02 | 0.02 | [ 0.98, 1.07] | 0.837 | [0.92, 1.08] | 0.992 | 1.002 | 2103 | 3042 |
| Daily pressure experienced | 0.90* | 0.04 | [ 0.82, 1.00] | 0.977 | [0.92, 1.08] | 0.302 | 1.000 | 4858 | 2932 |
| Daily pressure utilized (partner’s view) | 0.93 | 0.04 | [ 0.85, 1.02] | 0.940 | [0.92, 1.08] | 0.599 | 1.002 | 4006 | 2834 |
| Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Day | 0.96 | 0.06 | [ 0.85, 1.08] | 0.750 | [0.92, 1.08] | 0.700 | 1.002 | 8120 | 3138 |
| Own actionplan | 1.33*** | 0.06 | [ 1.22, 1.46] | 1.000 | [0.92, 1.08] | 0.000 | 1.000 | 5396 | 2975 |
| Partner actionplan | 1.08 | 0.04 | [ 0.99, 1.17] | 0.965 | [0.92, 1.08] | 0.538 | 1.001 | 6512 | 3309 |
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Conditional Between-Person Effects | |||||||||
| Mean persuasion experienced | 1.10 | 0.15 | [ 0.84, 1.45] | 0.758 | [0.92, 1.08] | 0.348 | 1.007 | 827 | 1540 |
| Mean persuasion utilized (partner’s view) | 1.08 | 0.15 | [ 0.82, 1.41] | 0.709 | [0.92, 1.08] | 0.378 | 1.007 | 833 | 1449 |
| Mean pressure experienced | 1.16 | 0.19 | [ 0.82, 1.59] | 0.808 | [0.92, 1.08] | 0.252 | 1.005 | 1105 | 2109 |
| Mean pressure utilized (partner’s view) | 0.96 | 0.16 | [ 0.69, 1.34] | 0.599 | [0.92, 1.08] | 0.346 | 1.005 | 1101 | 2103 |
| Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Within-Person Effects | |||||||||
| Hu Daily persuasion experienced | 1.54*** | 0.11 | [ 1.36, 1.79] | 1.000 | [0.84, 1.20] | 0.000 | 1.000 | 2888 | 2680 |
| Hu Daily persuasion utilized (partner’s view) | 1.36*** | 0.08 | [ 1.21, 1.54] | 1.000 | [0.84, 1.20] | 0.013 | 1.002 | 4621 | 3272 |
| Hu Daily pressure experienced | 0.97 | 0.13 | [ 0.73, 1.28] | 0.592 | [0.84, 1.20] | 0.792 | 1.000 | 5005 | 2787 |
| Hu Daily pressure utilized (partner’s view) | 1.58** | 0.29 | [ 1.14, 2.42] | 0.997 | [0.84, 1.20] | 0.050 | 1.000 | 3681 | 1761 |
| Hu Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Day | 0.85 | 0.12 | [ 0.64, 1.14] | 0.862 | [0.84, 1.20] | 0.537 | 1.001 | 7820 | 2875 |
| Hu Own actionplan | 8.95*** | 0.86 | [ 7.45, 10.88] | 1.000 | [0.84, 1.20] | 0.000 | 1.000 | 5652 | 3087 |
| Hu Partner actionplan | 1.22* | 0.12 | [ 1.01, 1.47] | 0.983 | [0.84, 1.20] | 0.420 | 1.003 | 5463 | 2654 |
| Hu Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Between-Person Effects | |||||||||
| Hu Mean persuasion experienced | 1.51 | 0.48 | [ 0.81, 2.77] | 0.896 | [0.84, 1.20] | 0.196 | 1.003 | 1235 | 2156 |
| Hu Mean persuasion utilized (partner’s view) | 1.55 | 0.48 | [ 0.84, 2.82] | 0.918 | [0.84, 1.20] | 0.180 | 1.003 | 1187 | 1565 |
| Hu Mean pressure experienced | 0.35** | 0.13 | [ 0.17, 0.73] | 0.997 | [0.84, 1.20] | 0.010 | 1.004 | 1541 | 2601 |
| Hu Mean pressure utilized (partner’s view) | 0.58 | 0.21 | [ 0.28, 1.23] | 0.923 | [0.84, 1.20] | 0.127 | 1.001 | 1876 | 2284 |
| Hu Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| sd(Intercept) | 0.30 | 0.04 | [0.23, 0.40] | NA | NA | NA | 1.004 | 1351 | 2570 |
| sd(Hurdle Intercept) | 0.72 | 0.10 | [0.54, 0.96] | NA | NA | NA | 1.001 | 1349 | 2021 |
| sd(Daily persuasion experienced) | 0.12 | 0.02 | [0.08, 0.16] | NA | NA | NA | 1.001 | 1975 | 2626 |
| sd(Daily persuasion utilized (partner’s view)) | 0.08 | 0.02 | [0.04, 0.13] | NA | NA | NA | 1.001 | 2475 | 2465 |
| sd(Daily pressure experienced) | 0.07 | 0.06 | [0.00, 0.23] | NA | NA | NA | 1.003 | 1555 | 1609 |
| sd(Daily pressure utilized (partner’s view)) | 0.07 | 0.05 | [0.00, 0.19] | NA | NA | NA | 1.003 | 1878 | 1845 |
| sd(Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily persuasion experienced) | 0.22 | 0.08 | [0.07, 0.41] | NA | NA | NA | 1.003 | 1464 | 1606 |
| sd(Hu Daily persuasion utilized (partner’s view)) | 0.16 | 0.09 | [0.01, 0.34] | NA | NA | NA | 1.003 | 1315 | 1776 |
| sd(Hu Daily pressure experienced) | 0.17 | 0.16 | [0.01, 0.68] | NA | NA | NA | 1.000 | 2157 | 2330 |
| sd(Hu Daily pressure utilized (partner’s view)) | 0.23 | 0.21 | [0.01, 0.89] | NA | NA | NA | 1.001 | 1678 | 1995 |
| sd(Hu Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||
| sigma | 0.68 | 0.01 | [0.66, 0.70] | NA | NA | NA | 1.003 | 6538 | 2471 |
# Plot continuous part of model
variable <- c(
'(Intercept)',
'b_persuasion_self_cw',
'b_persuasion_partner_cw',
'b_pressure_self_cw',
'b_pressure_partner_cw'
)
plot(
bayestestR::p_direction(pa_sub, parameter = variable),
priors = TRUE
) + theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
plot(
bayestestR::rope(
pa_sub,
parameter = variable,
range = rope_range_continuous,
verbose = F,
ci = 1
)
) + theme_bw()
# Hurdle part of the model
variable <- c(
'b_hu_persuasion_self_cw',
'b_hu_persuasion_partner_cw',
'b_hu_pressure_self_cw',
'b_hu_pressure_partner_cw'
)
plot(
bayestestR::p_direction(pa_sub, parameter = variable),
priors = TRUE
) + theme_bw()
# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
verbose = FALSE
) + theme_bw()## Possible multicollinearity between b_persuasion_partner_cb and
## b_persuasion_self_cb (r = 0.78), b_pressure_self_cb and
## b_persuasion_self_cb (r = 0.79), b_pressure_partner_cb and
## b_persuasion_partner_cb (r = 0.77), b_hu_persuasion_partner_cb and
## b_hu_persuasion_self_cb (r = 0.71), b_hu_pressure_self_cb and
## b_hu_persuasion_self_cb (r = 0.76), b_hu_pressure_partner_cb and
## b_hu_persuasion_partner_cb (r = 0.75). This might lead to inappropriate
## results. See 'Details' in '?rope'.
Device Based MVPA
Lognormal Model
formula <- bf(
pa_obj ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
plan_self + plan_partner +
day + weartime_self_cw + weartime_self_cb +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "Intercept")
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix_nopushing))
)## Warning: Rows containing NAs were excluded from the model.# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summarize_brms(
pa_obj_log,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
rope_range = rope_range_log,
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack)## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 110.54*** | 6.14 | [99.67, 123.17] | 1.000 | [0.94, 1.07] | 0.000 | 1.002 | 744 | 1473 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | 1.03 | 0.01 | [ 1.00, 1.06] | 0.972 | [0.94, 1.07] | 0.992 | 1.002 | 2381 | 2780 |
| Daily persuasion utilized (partner’s view) | 1.02 | 0.02 | [ 0.99, 1.05] | 0.883 | [0.94, 1.07] | 0.999 | 1.000 | 3188 | 3051 |
| Daily pressure experienced | 0.95 | 0.03 | [ 0.89, 1.02] | 0.934 | [0.94, 1.07] | 0.654 | 1.000 | 4197 | 2748 |
| Daily pressure utilized (partner’s view) | 0.98 | 0.03 | [ 0.92, 1.05] | 0.690 | [0.94, 1.07] | 0.927 | 1.001 | 4373 | 2658 |
| Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Day | 0.97 | 0.03 | [ 0.91, 1.04] | 0.797 | [0.94, 1.07] | 0.838 | 1.001 | 5395 | 2627 |
| Own Actionplan | 1.07** | 0.03 | [ 1.02, 1.12] | 0.996 | [0.94, 1.07] | 0.498 | 1.001 | 4381 | 2969 |
| Partner Actionplan | 1.05* | 0.03 | [ 1.01, 1.11] | 0.988 | [0.94, 1.07] | 0.680 | 1.002 | 5206 | 3137 |
| Between-Person Effects | |||||||||
| Daily weartime | 1.00*** | 0.00 | [ 1.00, 1.00] | 1.000 | [0.94, 1.07] | 1.000 | 1.000 | 3903 | 2616 |
| Mean persuasion experienced | 1.09 | 0.13 | [ 0.85, 1.41] | 0.763 | [0.94, 1.07] | 0.304 | 1.003 | 765 | 959 |
| Mean persuasion utilized (partner’s view) | 1.02 | 0.13 | [ 0.79, 1.30] | 0.566 | [0.94, 1.07] | 0.382 | 1.003 | 751 | 1151 |
| Mean pressure experienced | 0.91 | 0.13 | [ 0.69, 1.22] | 0.738 | [0.94, 1.07] | 0.285 | 1.002 | 949 | 1514 |
| Mean pressure utilized (partner’s view) | 1.04 | 0.14 | [ 0.80, 1.37] | 0.618 | [0.94, 1.07] | 0.349 | 1.003 | 873 | 1354 |
| Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| Mean weartime | 1.00 | 0.00 | [ 1.00, 1.00] | 0.912 | [0.94, 1.07] | 1.000 | 1.002 | 4283 | 3057 |
| sd(Intercept) | 0.29 | 0.04 | [0.23, 0.38] | NA | NA | NA | 1.001 | 908 | 1500 |
| sd(Daily persuasion experienced) | 0.05 | 0.01 | [0.02, 0.08] | NA | NA | NA | 1.000 | 1601 | 1315 |
| sd(Daily persuasion utilized (partner’s view)) | 0.05 | 0.01 | [0.02, 0.09] | NA | NA | NA | 1.001 | 1940 | 2132 |
| sd(Daily pressure experienced) | 0.04 | 0.04 | [0.00, 0.14] | NA | NA | NA | 1.004 | 1369 | 1870 |
| sd(Daily pressure utilized (partner’s view)) | 0.04 | 0.03 | [0.00, 0.11] | NA | NA | NA | 1.000 | 2374 | 2184 |
| sd(Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||
| sd(Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sigma | 0.58 | 0.01 | [0.56, 0.59] | NA | NA | NA | 1.000 | 5874 | 2881 |
plot(
bayestestR::p_direction(pa_obj_log),
priors = TRUE
) +
coord_cartesian(xlim = c(-3, 3)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Possible multicollinearity between b_persuasion_partner_cb and
## b_persuasion_self_cb (r = 0.88), b_pressure_self_cb and
## b_persuasion_self_cb (r = 0.84), b_pressure_partner_cb and
## b_persuasion_self_cb (r = 0.78), b_pressure_self_cb and
## b_persuasion_partner_cb (r = 0.76), b_pressure_partner_cb and
## b_persuasion_partner_cb (r = 0.87). This might lead to inappropriate
## results. See 'Details' in '?rope'.
Affect
Gaussian
formula <- bf(
aff ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
plan_self + plan_partner +
day +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss", suffix_nopushing))
)## Warning: Rows containing NAs were excluded from the model.summarize_brms(
mood_gauss,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = F) %>%
print_df(rows_to_pack = rows_to_pack)| Est. | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 3.67*** | 0.11 | [ 3.46, 3.89] | 1.000 | [-0.11, 0.11] | 0.000 | 1.009 | 419 | 981 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | 0.00 | 0.02 | [-0.04, 0.05] | 0.560 | [-0.11, 0.11] | 1.000 | 1.001 | 3858 | 3241 |
| Daily persuasion utilized (partner’s view) | 0.03 | 0.02 | [-0.01, 0.08] | 0.909 | [-0.11, 0.11] | 0.997 | 1.000 | 2990 | 2740 |
| Daily pressure experienced | -0.03 | 0.05 | [-0.14, 0.07] | 0.725 | [-0.11, 0.11] | 0.943 | 1.003 | 4499 | 2684 |
| Daily pressure utilized (partner’s view) | 0.00 | 0.05 | [-0.11, 0.10] | 0.522 | [-0.11, 0.11] | 0.967 | 1.001 | 4772 | 3218 |
| Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Day | 0.25*** | 0.06 | [ 0.14, 0.36] | 1.000 | [-0.11, 0.11] | 0.007 | 1.000 | 7438 | 2668 |
| Own Actionplan | 0.11** | 0.04 | [ 0.04, 0.18] | 0.998 | [-0.11, 0.11] | 0.569 | 1.001 | 5645 | 3193 |
| Partner Actionplan | -0.01 | 0.04 | [-0.08, 0.07] | 0.584 | [-0.11, 0.11] | 0.997 | 1.000 | 6635 | 3351 |
| Between-Person Effects | |||||||||
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion experienced | 0.40 | 0.27 | [-0.10, 0.94] | 0.943 | [-0.11, 0.11] | 0.099 | 1.007 | 359 | 598 |
| Mean persuasion utilized (partner’s view) | 0.33 | 0.27 | [-0.19, 0.85] | 0.897 | [-0.11, 0.11] | 0.160 | 1.007 | 362 | 497 |
| Mean pressure experienced | -0.34 | 0.29 | [-0.91, 0.21] | 0.887 | [-0.11, 0.11] | 0.158 | 1.007 | 430 | 640 |
| Mean pressure utilized (partner’s view) | -0.28 | 0.29 | [-0.84, 0.28] | 0.840 | [-0.11, 0.11] | 0.206 | 1.007 | 407 | 616 |
| Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Intercept) | 0.59 | 0.08 | [0.47, 0.78] | NA | NA | NA | 1.004 | 602 | 1150 |
| sd(Daily persuasion experienced) | 0.05 | 0.03 | [0.00, 0.11] | NA | NA | NA | 1.003 | 982 | 1807 |
| sd(Daily persuasion utilized (partner’s view)) | 0.08 | 0.03 | [0.02, 0.13] | NA | NA | NA | 1.007 | 1022 | 681 |
| sd(Daily pressure experienced) | 0.06 | 0.06 | [0.00, 0.24] | NA | NA | NA | 1.002 | 1830 | 2395 |
| sd(Daily pressure utilized (partner’s view)) | 0.07 | 0.06 | [0.00, 0.24] | NA | NA | NA | 1.000 | 1764 | 1629 |
| sd(Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||
| sd(Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sigma | 0.96 | 0.01 | [0.94, 0.98] | NA | NA | NA | 1.000 | 7952 | 2848 |
plot(
bayestestR::p_direction(mood_gauss),
priors = TRUE
) +
coord_cartesian(xlim = c(-3, 3)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Possible multicollinearity between b_pressure_partner_cb and
## b_persuasion_self_cb (r = 0.87), b_pressure_self_cb and
## b_persuasion_partner_cb (r = 0.87), b_pressure_partner_cb and
## b_pressure_self_cb (r = 0.82). This might lead to inappropriate results.
## See 'Details' in '?rope'.
Reactance
Ordinal
df_double$reactance_ordinal <- factor(df_double$reactance,
levels = 0:5,
ordered = TRUE)
formula <- bf(
reactance_ordinal ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
plan_self + plan_partner +
day +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = cumulative() # HURDLE_CUMULATIVE
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
reactance_ordinal <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::cumulative(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777
, file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix_nopushing))
)## Warning: Rows containing NAs were excluded from the model.summarize_brms(
reactance_ordinal,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
rope_range = c(-0.18, 0.18),
model_rows_fixed = model_rows_fixed_ordinal,
model_rows_random = model_rows_random_ordinal,
model_rownames_fixed = model_rownames_fixed_ordinal,
model_rownames_random = model_rownames_random_ordinal,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack_ordinal)| OR | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercepts | |||||||||
| Intercept | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Intercept[1] | 3.46*** | 0.99 | [ 1.99, 6.09] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 4277 | 3245 |
| Intercept[2] | 7.31*** | 2.21 | [ 4.11, 13.20] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 4239 | 3268 |
| Intercept[3] | 19.71*** | 6.39 | [ 10.86, 37.64] | 1.000 | [0.84, 1.20] | 0.000 | 1.000 | 4388 | 3235 |
| Intercept[4] | 85.54*** | 32.07 | [ 40.85, 177.13] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 4946 | 3129 |
| Intercept[5] | 2773.94*** | 1761.22 | [844.52, 11920.75] | 1.000 | [0.84, 1.20] | 0.000 | 1.002 | 7124 | 3328 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | 0.85* | 0.07 | [ 0.72, 0.99] | 0.980 | [0.84, 1.20] | 0.580 | 1.001 | 4159 | 3362 |
| Daily persuasion utilized (partner’s view) | 1.01 | 0.09 | [ 0.85, 1.21] | 0.564 | [0.84, 1.20] | 0.948 | 1.000 | 3583 | 2527 |
| Daily pressure experienced | 2.02** | 0.37 | [ 1.34, 2.84] | 0.999 | [0.84, 1.20] | 0.009 | 1.003 | 2038 | 2328 |
| Daily pressure utilized (partner’s view) | 1.13 | 0.24 | [ 0.69, 1.77] | 0.718 | [0.84, 1.20] | 0.527 | 1.002 | 2561 | 2272 |
| Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Day | 1.64 | 0.53 | [ 0.88, 3.08] | 0.934 | [0.84, 1.20] | 0.148 | 1.000 | 5394 | 2906 |
| Own Actionplan | 1.12 | 0.26 | [ 0.71, 1.77] | 0.688 | [0.84, 1.20] | 0.504 | 1.001 | 4935 | 3695 |
| Partner Actionplan | 0.86 | 0.19 | [ 0.55, 1.36] | 0.745 | [0.84, 1.20] | 0.476 | 1.002 | 5147 | 3304 |
| Between-Person Effects | |||||||||
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion experienced | 1.11 | 0.53 | [ 0.44, 2.92] | 0.582 | [0.84, 1.20] | 0.278 | 1.000 | 1543 | 2308 |
| Mean persuasion utilized (partner’s view) | 0.80 | 0.42 | [ 0.29, 2.19] | 0.673 | [0.84, 1.20] | 0.251 | 1.002 | 1841 | 2306 |
| Mean pressure experienced | 3.27* | 1.69 | [ 1.20, 9.52] | 0.991 | [0.84, 1.20] | 0.020 | 1.000 | 1663 | 2685 |
| Mean pressure utilized (partner’s view) | 1.21 | 0.69 | [ 0.41, 3.57] | 0.627 | [0.84, 1.20] | 0.232 | 1.001 | 1989 | 2288 |
| Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Intercept) | 0.78 | 0.19 | [0.46, 1.26] | NA | NA | NA | 1.000 | 1675 | 2397 |
| sd(Daily persuasion experienced) | 0.17 | 0.13 | [0.01, 0.44] | NA | NA | NA | 1.005 | 785 | 1414 |
| sd(Daily persuasion utilized (partner’s view)) | 0.19 | 0.12 | [0.01, 0.47] | NA | NA | NA | 1.001 | 1330 | 1750 |
| sd(Daily pressure experienced) | 0.47 | 0.23 | [0.07, 1.01] | NA | NA | NA | 1.005 | 985 | 890 |
| sd(Daily pressure utilized (partner’s view)) | 0.30 | 0.29 | [0.01, 1.37] | NA | NA | NA | 1.001 | 1106 | 2024 |
| sd(Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||
| sd(Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| disc | 1.00 | 0.00 | [1.00, 1.00] | NA | NA | NA | NA | NA | NA |
plot(
bayestestR::p_direction(reactance_ordinal),
priors = TRUE
) +
coord_cartesian(xlim = c(-6, 6)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
## = 0.79), b_Intercept[4] and b_Intercept[3] (r = 0.86),
## b_pressure_self_cb and b_persuasion_self_cb (r = 0.79),
## b_pressure_partner_cb and b_persuasion_partner_cb (r = 0.84). This might
## lead to inappropriate results. See 'Details' in '?rope'.
conditional_spaghetti(
reactance_ordinal,
effects = c('persuasion_self_cw', 'pressure_self_cw')
, group_var = 'coupleID'
#, n_groups = 15
, plot_full_range = T
)\(persuasion_self_cw
<img
src="03_SensitivityPushingSeparate_files/figure-html/report_reactance_ordinal-3.png"
width="2400" />\)pressure_self_cw
Binary
introduce_binary_reactance <- function(data) {
data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
return(data)
}
df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
for (i in seq_along(implist)) {
implist[[i]] <- introduce_binary_reactance(implist[[i]])
}
}
formula <- bf(
is_reactance ~
persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw +
persuasion_self_cb + persuasion_partner_cb +
pressure_self_cb + pressure_partner_cb +
plan_self + plan_partner +
day +
# Random effects
(persuasion_self_cw + persuasion_partner_cw +
pressure_self_cw + pressure_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5)
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance", suffix_nopushing))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.summarize_brms(
is_reactance,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack)## Warning: There were 4 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup| OR | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 0.31*** | 0.10 | [0.17, 0.57] | 1.000 | [0.83, 1.20] | 0.000 | 1.001 | 3635 | 3242 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | 0.85 | 0.08 | [0.70, 1.02] | 0.962 | [0.83, 1.20] | 0.605 | 1.002 | 3818 | 2660 |
| Daily persuasion utilized (partner’s view) | 1.11 | 0.14 | [0.87, 1.49] | 0.797 | [0.83, 1.20] | 0.706 | 1.000 | 3184 | 2493 |
| Daily pressure experienced | 2.06* | 0.53 | [1.20, 4.07] | 0.992 | [0.83, 1.20] | 0.021 | 1.000 | 2896 | 2535 |
| Daily pressure utilized (partner’s view) | 1.31 | 0.48 | [0.62, 3.89] | 0.775 | [0.83, 1.20] | 0.290 | 1.003 | 1265 | 560 |
| Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Day | 1.87 | 0.68 | [0.91, 3.78] | 0.957 | [0.83, 1.20] | 0.094 | 1.003 | 5369 | 3275 |
| Own Actionplan | 1.18 | 0.33 | [0.66, 2.05] | 0.718 | [0.83, 1.20] | 0.418 | 1.001 | 4143 | 1857 |
| Partner Actionplan | 0.87 | 0.23 | [0.52, 1.47] | 0.696 | [0.83, 1.20] | 0.450 | 1.000 | 4581 | 3096 |
| Between-Person Effects | |||||||||
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion experienced | 1.89 | 1.10 | [0.63, 6.01] | 0.870 | [0.83, 1.20] | 0.140 | 1.002 | 1752 | 2394 |
| Mean persuasion utilized (partner’s view) | 1.04 | 0.63 | [0.33, 3.68] | 0.524 | [0.83, 1.20] | 0.243 | 1.001 | 1993 | 2757 |
| Mean pressure experienced | 18.30*** | 16.32 | [3.32, 127.88] | 1.000 | [0.83, 1.20] | 0.001 | 1.000 | 2763 | 2477 |
| Mean pressure utilized (partner’s view) | 1.13 | 1.09 | [0.17, 7.51] | 0.555 | [0.83, 1.20] | 0.152 | 1.001 | 2116 | 2956 |
| Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Intercept) | 1.21 | 0.25 | [0.80, 1.79] | NA | NA | NA | 1.003 | 1504 | 1765 |
| sd(Daily persuasion experienced) | 0.19 | 0.13 | [0.01, 0.49] | NA | NA | NA | 1.007 | 987 | 1080 |
| sd(Daily persuasion utilized (partner’s view)) | 0.44 | 0.17 | [0.13, 0.87] | NA | NA | NA | 1.005 | 1425 | 1011 |
| sd(Daily pressure experienced) | 0.74 | 0.51 | [0.04, 2.04] | NA | NA | NA | 1.002 | 839 | 1539 |
| sd(Daily pressure utilized (partner’s view)) | 0.72 | 0.63 | [0.05, 2.65] | NA | NA | NA | 1.002 | 1040 | 969 |
| sd(Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||
| sd(Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA |
plot(
bayestestR::p_direction(is_reactance),
priors = TRUE
) +
coord_cartesian(xlim = c(-6, 6)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
Report All Models
process_model_component <- function(obj) {
# Convert to string, modify, and evaluate
name <- deparse(substitute(obj))
if (report_hurdle) name <- paste0(name, '_hu')
if (report_ordinal) name <- paste0(name, '_ordinal')
return(get(name, envir = parent.frame()))
}
all_models <- report_side_by_side(
pa_sub,
pa_obj_log,
mood_gauss,
reactance_ordinal,
is_reactance,
stats_to_report = c('CI', 'pd'),
model_rows_random = process_model_component(model_rows_random),
model_rows_fixed = process_model_component(model_rows_fixed),
model_rownames_random = process_model_component(model_rownames_random),
model_rownames_fixed = process_model_component(model_rownames_fixed)
) [1] “pa_sub”
## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.[1] “pa_obj_log”
## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”
## Warning: There were 4 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup# pretty printing
summary_all_models <- all_models %>%
print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
add_header_above(
c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),
"Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5),
"Mood Gaussian" = (length(all_models) / 5),
"Reactance Ordinal" = (length(all_models) / 5),
"Reactance Dichotome" = (length(all_models) / 5))
)
export_xlsx(
summary_all_models,
rows_to_pack = process_model_component(rows_to_pack),
file.path("Output", paste0("AllModels", suffix_nopushing, ".xlsx")),
merge_option = 'header',
simplify_2nd_row = TRUE,
line_above_rows = c(1,2),
line_below_rows = c(-1)
)##
## Attaching package: 'rvest'## The following object is masked from 'package:readr':
##
## guess_encoding| exp(Est.) pa_sub | 95% CI pa_sub | pd pa_sub | exp(Est.) pa_obj_log | 95% CI pa_obj_log | pd pa_obj_log | Est. mood_gauss | 95% CI mood_gauss | pd mood_gauss | OR reactance_ordinal | 95% CI reactance_ordinal | pd reactance_ordinal | OR is_reactance | 95% CI is_reactance | pd is_reactance | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Intercept | 38.46*** | [33.40, 44.04] | 1.000 | 110.54*** | [99.67, 123.17] | 1.000 | 3.67*** | [ 3.46, 3.89] | 1.000 | NA | NA | NA | 0.31*** | [0.17, 0.57] | 1.000 |
| Hurdle Intercept | 0.25*** | [ 0.19, 0.35] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Conditional Within-Person Effects | |||||||||||||||
| Daily persuasion experienced | 1.02 | [ 0.97, 1.08] | 0.817 | 1.03 | [ 1.00, 1.06] | 0.972 | 0.00 | [-0.04, 0.05] | 0.560 | 0.85* | [ 0.72, 0.99] | 0.980 | 0.85 | [0.70, 1.02] | 0.962 |
| Daily persuasion utilized (partner’s view) | 1.02 | [ 0.98, 1.07] | 0.837 | 1.02 | [ 0.99, 1.05] | 0.883 | 0.03 | [-0.01, 0.08] | 0.909 | 1.01 | [ 0.85, 1.21] | 0.564 | 1.11 | [0.87, 1.49] | 0.797 |
| Daily pressure experienced | 0.90* | [ 0.82, 1.00] | 0.977 | 0.95 | [ 0.89, 1.02] | 0.934 | -0.03 | [-0.14, 0.07] | 0.725 | 2.02** | [ 1.34, 2.84] | 0.999 | 2.06* | [1.20, 4.07] | 0.992 |
| Daily pressure utilized (partner’s view) | 0.93 | [ 0.85, 1.02] | 0.940 | 0.98 | [ 0.92, 1.05] | 0.690 | 0.00 | [-0.11, 0.10] | 0.522 | 1.13 | [ 0.69, 1.77] | 0.718 | 1.31 | [0.62, 3.89] | 0.775 |
| Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Day | 0.96 | [ 0.85, 1.08] | 0.750 | 0.97 | [ 0.91, 1.04] | 0.797 | 0.25*** | [ 0.14, 0.36] | 1.000 | 1.64 | [ 0.88, 3.08] | 0.934 | 1.87 | [0.91, 3.78] | 0.957 |
| Own actionplan | 1.33*** | [ 1.22, 1.46] | 1.000 | 1.07** | [ 1.02, 1.12] | 0.996 | 0.11** | [ 0.04, 0.18] | 0.998 | 1.12 | [ 0.71, 1.77] | 0.688 | 1.18 | [0.66, 2.05] | 0.718 |
| Partner actionplan | 1.08 | [ 0.99, 1.17] | 0.965 | 1.05* | [ 1.01, 1.11] | 0.988 | -0.01 | [-0.08, 0.07] | 0.584 | 0.86 | [ 0.55, 1.36] | 0.745 | 0.87 | [0.52, 1.47] | 0.696 |
| Daily weartime | NA | NA | NA | 1.00*** | [ 1.00, 1.00] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Conditional Between-Person Effects | |||||||||||||||
| Mean persuasion experienced | 1.10 | [ 0.84, 1.45] | 0.758 | 1.09 | [ 0.85, 1.41] | 0.763 | 0.40 | [-0.10, 0.94] | 0.943 | 1.11 | [ 0.44, 2.92] | 0.582 | 1.89 | [0.63, 6.01] | 0.870 |
| Mean persuasion utilized (partner’s view) | 1.08 | [ 0.82, 1.41] | 0.709 | 1.02 | [ 0.79, 1.30] | 0.566 | 0.33 | [-0.19, 0.85] | 0.897 | 0.80 | [ 0.29, 2.19] | 0.673 | 1.04 | [0.33, 3.68] | 0.524 |
| Mean pressure experienced | 1.16 | [ 0.82, 1.59] | 0.808 | 0.91 | [ 0.69, 1.22] | 0.738 | -0.34 | [-0.91, 0.21] | 0.887 | 3.27* | [ 1.20, 9.52] | 0.991 | 18.30*** | [3.32, 127.88] | 1.000 |
| Mean pressure utilized (partner’s view) | 0.96 | [ 0.69, 1.34] | 0.599 | 1.04 | [ 0.80, 1.37] | 0.618 | -0.28 | [-0.84, 0.28] | 0.840 | 1.21 | [ 0.41, 3.57] | 0.627 | 1.13 | [0.17, 7.51] | 0.555 |
| Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean weartime | NA | NA | NA | 1.00 | [ 1.00, 1.00] | 0.912 | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Within-Person Effects | |||||||||||||||
| Hu Daily persuasion experienced | 1.54*** | [ 1.36, 1.79] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily persuasion utilized (partner’s view) | 1.36*** | [ 1.21, 1.54] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pressure experienced | 0.97 | [ 0.73, 1.28] | 0.592 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pressure utilized (partner’s view) | 1.58** | [ 1.14, 2.42] | 0.997 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Day | 0.85 | [ 0.64, 1.14] | 0.862 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Own actionplan | 8.95*** | [ 7.45, 10.88] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Partner actionplan | 1.22* | [ 1.01, 1.47] | 0.983 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Between-Person Effects | |||||||||||||||
| Hu Mean persuasion experienced | 1.51 | [ 0.81, 2.77] | 0.896 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean persuasion utilized (partner’s view) | 1.55 | [ 0.84, 2.82] | 0.918 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pressure experienced | 0.35** | [ 0.17, 0.73] | 0.997 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pressure utilized (partner’s view) | 0.58 | [ 0.28, 1.23] | 0.923 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pushing experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pushing utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||||||||
| sd(Intercept) | 0.30 | [0.23, 0.40] | NA | 0.29 | [0.23, 0.38] | NA | 0.59 | [0.47, 0.78] | NA | 0.78 | [0.46, 1.26] | NA | 1.21 | [0.80, 1.79] | NA |
| sd(Hurdle Intercept) | 0.72 | [0.54, 0.96] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion experienced) | 0.12 | [0.08, 0.16] | NA | 0.05 | [0.02, 0.08] | NA | 0.05 | [0.00, 0.11] | NA | 0.17 | [0.01, 0.44] | NA | 0.19 | [0.01, 0.49] | NA |
| sd(Daily persuasion utilized (partner’s view)) | 0.08 | [0.04, 0.13] | NA | 0.05 | [0.02, 0.09] | NA | 0.08 | [0.02, 0.13] | NA | 0.19 | [0.01, 0.47] | NA | 0.44 | [0.13, 0.87] | NA |
| sd(Daily pressure experienced) | 0.07 | [0.00, 0.23] | NA | 0.04 | [0.00, 0.14] | NA | 0.06 | [0.00, 0.24] | NA | 0.47 | [0.07, 1.01] | NA | 0.74 | [0.04, 2.04] | NA |
| sd(Daily pressure utilized (partner’s view)) | 0.07 | [0.00, 0.19] | NA | 0.04 | [0.00, 0.11] | NA | 0.07 | [0.00, 0.24] | NA | 0.30 | [0.01, 1.37] | NA | 0.72 | [0.05, 2.65] | NA |
| sd(Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily persuasion experienced) | 0.22 | [0.07, 0.41] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily persuasion utilized (partner’s view)) | 0.16 | [0.01, 0.34] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pressure experienced) | 0.17 | [0.01, 0.68] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pressure utilized (partner’s view)) | 0.23 | [0.01, 0.89] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pushing experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pushing utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||||||||
| sigma | 0.68 | [0.66, 0.70] | NA | 0.58 | [0.56, 0.59] | NA | 0.96 | [0.94, 0.98] | NA | NA | NA | NA | NA | NA | NA |
ONLY PUSHING
Hurdle Lognormal Model
formula <- bf(
pa_sub ~
pushing_self_cw + pushing_partner_cw +
pushing_self_cb + pushing_partner_cb +
#plan_self + plan_partner +
day +
# Random effects
(1 + pushing_self_cw + pushing_partner_cw | coupleID),
hu = ~ pushing_self_cw + pushing_partner_cw +
pushing_self_cb + pushing_partner_cb +
#plan_self + plan_partner +
day +
# Random effects
(1 + pushing_self_cw + pushing_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
, brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
, brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = hurdle_lognormal()
#)
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_sub <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::hurdle_lognormal(),
#family = brms::hurdle_negbinomial(),
#family = brms::hurdle_poisson(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 42,
file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix_onlypushing))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
pa_sub,
stats_to_report = c('pd','ROPE'),
rope_range = c(-0.18, 0.18),
model_rows_fixed = model_rows_fixed_hu,
model_rows_random = model_rows_random_hu,
model_rownames_fixed = model_rownames_fixed_hu,
model_rownames_random = model_rownames_random_hu,
exponentiate = T) ## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)
summary_pa_sub_continuous <- summarize_brms(
pa_sub,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
rope_range = rope_range_continuous,
model_rows_fixed = model_rows_fixed_hu,
model_rows_random = model_rows_random_hu,
model_rownames_fixed = model_rownames_fixed_hu,
model_rownames_random = model_rownames_random_hu,
exponentiate = T) ## Warning in summarize_brms(pa_sub, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous
columns_to_replace <- c("ROPE", "inside ROPE")
summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <-
summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]
# Print the updated dataframe
summary_pa_sub %>%
print_df(rows_to_pack = rows_to_pack_hu)| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 57.28*** | 4.53 | [49.10, 66.97] | 1.000 | [0.93, 1.08] | 0.000 | 1.002 | 1327 | 2040 |
| Hurdle Intercept | 3.79*** | 0.98 | [ 2.27, 6.26] | 1.000 | [0.84, 1.20] | 0.000 | 1.005 | 1088 | 2463 |
| Conditional Within-Person Effects | |||||||||
| Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing experienced | 0.97 | 0.03 | [ 0.91, 1.03] | 0.856 | [0.93, 1.08] | 0.931 | 1.002 | 4427 | 2662 |
| Daily pushing utilized (partner’s view) | 0.97 | 0.03 | [ 0.90, 1.03] | 0.847 | [0.93, 1.08] | 0.907 | 1.001 | 3679 | 2873 |
| Day | 0.91 | 0.07 | [ 0.79, 1.06] | 0.894 | [0.93, 1.08] | 0.398 | 1.000 | 8744 | 2600 |
| Own actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Partner actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Conditional Between-Person Effects | |||||||||
| Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing experienced | 1.32* | 0.18 | [ 1.00, 1.71] | 0.977 | [0.93, 1.08] | 0.060 | 1.005 | 1217 | 2176 |
| Mean pushing utilized (partner’s view) | 1.18 | 0.15 | [ 0.90, 1.53] | 0.883 | [0.93, 1.08] | 0.222 | 1.005 | 1265 | 2076 |
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Within-Person Effects | |||||||||
| Hu Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pushing experienced | 1.55** | 0.23 | [ 1.17, 2.24] | 0.998 | [0.84, 1.20] | 0.035 | 1.002 | 3127 | 2531 |
| Hu Daily pushing utilized (partner’s view) | 1.84*** | 0.30 | [ 1.38, 2.68] | 1.000 | [0.84, 1.20] | 0.002 | 1.001 | 2766 | 2989 |
| Hu Day | 0.61 | 0.15 | [ 0.37, 1.00] | 0.974 | [0.84, 1.20] | 0.100 | 1.003 | 7632 | 2768 |
| Hu Own actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Partner actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Between-Person Effects | |||||||||
| Hu Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pushing experienced | 0.46** | 0.14 | [ 0.26, 0.81] | 0.997 | [0.84, 1.20] | 0.018 | 1.001 | 2675 | 2824 |
| Hu Mean pushing utilized (partner’s view) | 0.62 | 0.18 | [ 0.33, 1.08] | 0.952 | [0.84, 1.20] | 0.136 | 1.001 | 2568 | 2699 |
| Hu Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||
| sd(Intercept) | 0.37 | 0.06 | [0.29, 0.50] | NA | NA | NA | 1.002 | 1283 | 2370 |
| sd(Hurdle Intercept) | 1.25 | 0.19 | [0.93, 1.74] | NA | NA | NA | 1.002 | 1448 | 2426 |
| sd(Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing experienced) | 0.06 | 0.04 | [0.00, 0.16] | NA | NA | NA | 1.001 | 870 | 1486 |
| sd(Daily pushing utilized (partner’s view)) | 0.09 | 0.04 | [0.01, 0.18] | NA | NA | NA | 1.004 | 977 | 1142 |
| sd(Hu Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pushing experienced) | 0.46 | 0.24 | [0.07, 1.02] | NA | NA | NA | 1.000 | 947 | 1689 |
| sd(Hu Daily pushing utilized (partner’s view)) | 0.51 | 0.24 | [0.10, 1.13] | NA | NA | NA | 1.001 | 1088 | 1582 |
| Additional Parameters | |||||||||
| sigma | 0.67 | 0.02 | [0.64, 0.71] | NA | NA | NA | 1.002 | 6374 | 2773 |
# Plot continuous part of model
variable <- c(
'(Intercept)',
'b_pushing_self_cw',
'b_pushing_partner_cw'
)
plot(
bayestestR::p_direction(pa_sub, parameter = variable),
priors = TRUE
) + theme_bw()
plot(
bayestestR::rope(
pa_sub,
parameter = variable,
range = rope_range_continuous,
verbose = F,
ci = 1
)
) + theme_bw()
# Hurdle part of the model
variable <- c(
'b_hu_pushing_self_cw',
'b_hu_pushing_partner_cw'
)
plot(
bayestestR::p_direction(pa_sub, parameter = variable),
priors = TRUE
) + theme_bw()
# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
verbose = FALSE
) + theme_bw()## Possible multicollinearity between b_pushing_partner_cb and
## b_pushing_self_cb (r = 0.77). This might lead to inappropriate results.
## See 'Details' in '?rope'.
Device Based MVPA
Lognormal Model
formula <- bf(
pa_obj ~
pushing_self_cw + pushing_partner_cw +
pushing_self_cb + pushing_partner_cb +
#plan_self + plan_partner +
day + weartime_self_cw + weartime_self_cb +
# Random effects
(pushing_self_cw + pushing_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 50)", class = "Intercept")
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = lognormal()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
pa_obj_log <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = lognormal(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix_onlypushing))
)## Warning: Rows containing NAs were excluded from the model.# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)
summarize_brms(
pa_obj_log,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
rope_range = rope_range_log,
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack)## Warning in summarize_brms(pa_obj_log, stats_to_report = c("CI", "SE", "pd", :
## Coefficients were exponentiated. Double check if this was intended.| exp(Est.) | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 124.24*** | 7.52 | [109.69, 139.58] | 1.000 | [0.94, 1.06] | 0.000 | 1.002 | 967 | 1634 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing experienced | 1.03 | 0.03 | [ 0.97, 1.08] | 0.807 | [0.94, 1.06] | 0.914 | 1.001 | 3221 | 2733 |
| Daily pushing utilized (partner’s view) | 1.02 | 0.02 | [ 0.97, 1.06] | 0.797 | [0.94, 1.06] | 0.979 | 1.000 | 5065 | 2748 |
| Day | 0.97 | 0.05 | [ 0.87, 1.08] | 0.718 | [0.94, 1.06] | 0.654 | 1.000 | 7350 | 3336 |
| Own Actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Partner Actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Between-Person Effects | |||||||||
| Daily weartime | 1.00** | 0.00 | [ 1.00, 1.00] | 0.998 | [0.94, 1.06] | 1.000 | 1.001 | 4680 | 2898 |
| Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing experienced | 0.97 | 0.06 | [ 0.86, 1.10] | 0.680 | [0.94, 1.06] | 0.635 | 1.001 | 2223 | 2645 |
| Mean pushing utilized (partner’s view) | 1.04 | 0.07 | [ 0.92, 1.18] | 0.739 | [0.94, 1.06] | 0.578 | 1.002 | 2058 | 2692 |
| Random Effects | |||||||||
| Mean weartime | 1.00 | 0.00 | [ 1.00, 1.00] | 0.812 | [0.94, 1.06] | 1.000 | 1.000 | 3879 | 3246 |
| sd(Intercept) | 0.31 | 0.04 | [0.24, 0.42] | NA | NA | NA | 1.005 | 890 | 1790 |
| sd(Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing experienced) | 0.09 | 0.04 | [0.02, 0.17] | NA | NA | NA | 1.004 | 887 | 821 |
| Additional Parameters | |||||||||
| sd(Daily pushing utilized (partner’s view)) | 0.03 | 0.03 | [0.00, 0.10] | NA | NA | NA | 1.001 | 1398 | 1731 |
| sigma | 0.54 | 0.01 | [0.52, 0.57] | NA | NA | NA | 1.001 | 6040 | 3037 |
plot(
bayestestR::p_direction(pa_obj_log),
priors = TRUE
) +
coord_cartesian(xlim = c(-3, 3)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
Affect
Gaussian
formula <- bf(
aff ~
pushing_self_cw + pushing_partner_cw +
pushing_self_cb + pushing_partner_cb +
#plan_self + plan_partner +
day +
# Random effects
(pushing_self_cw + pushing_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 5)", class = "b")
,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
, brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = gaussian()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
mood_gauss <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = gaussian(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("mood_gauss", suffix_onlypushing))
)## Warning: Rows containing NAs were excluded from the model.summarize_brms(
mood_gauss,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = F) %>%
print_df(rows_to_pack = rows_to_pack)| Est. | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 3.78*** | 0.12 | [ 3.56, 4.00] | 1.000 | [-0.11, 0.11] | 0.000 | 1.013 | 560 | 1111 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing experienced | 0.02 | 0.04 | [-0.05, 0.08] | 0.670 | [-0.11, 0.11] | 0.995 | 1.002 | 2823 | 2591 |
| Daily pushing utilized (partner’s view) | 0.09* | 0.04 | [ 0.01, 0.17] | 0.984 | [-0.11, 0.11] | 0.715 | 1.002 | 1597 | 2207 |
| Day | 0.28*** | 0.09 | [ 0.11, 0.46] | 1.000 | [-0.11, 0.11] | 0.023 | 1.001 | 3871 | 3088 |
| Own Actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Partner Actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Between-Person Effects | |||||||||
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing experienced | 0.12 | 0.11 | [-0.08, 0.32] | 0.880 | [-0.11, 0.11] | 0.433 | 1.003 | 1269 | 1954 |
| Mean pushing utilized (partner’s view) | 0.08 | 0.11 | [-0.12, 0.30] | 0.786 | [-0.11, 0.11] | 0.554 | 1.002 | 1216 | 1958 |
| Random Effects | |||||||||
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Intercept) | 0.61 | 0.08 | [0.48, 0.79] | NA | NA | NA | 1.011 | 648 | 1448 |
| sd(Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing experienced) | 0.07 | 0.04 | [0.00, 0.16] | NA | NA | NA | 1.000 | 1409 | 1126 |
| Additional Parameters | |||||||||
| sd(Daily pushing utilized (partner’s view)) | 0.14 | 0.04 | [0.06, 0.24] | NA | NA | NA | 1.002 | 1329 | 828 |
| sigma | 0.91 | 0.02 | [0.87, 0.94] | NA | NA | NA | 1.001 | 3887 | 2536 |
plot(
bayestestR::p_direction(mood_gauss),
priors = TRUE
) +
coord_cartesian(xlim = c(-3, 3)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
conditional_spaghetti(
mood_gauss,
effects = c('pushing_partner_cw'),
group_var = 'coupleID',
plot_full_range = TRUE
)$pushing_partner_cw
Reactance
Ordinal
df_double$reactance_ordinal <- factor(df_double$reactance,
levels = 0:5,
ordered = TRUE)
formula <- bf(
reactance_ordinal ~
pushing_self_cw + pushing_partner_cw +
pushing_self_cb + pushing_partner_cb +
#plan_self + plan_partner +
day +
# Random effects
(pushing_self_cw + pushing_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = cumulative() # HURDLE_CUMULATIVE
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
reactance_ordinal <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::cumulative(),
#control = list(adapt_delta = 0.95),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777
, file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix_onlypushing))
)## Warning: Rows containing NAs were excluded from the model.summarize_brms(
reactance_ordinal,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
rope_range = c(-0.18, 0.18),
model_rows_fixed = model_rows_fixed_ordinal,
model_rows_random = model_rows_random_ordinal,
model_rownames_fixed = model_rownames_fixed_ordinal,
model_rownames_random = model_rownames_random_ordinal,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack_ordinal)| OR | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercepts | |||||||||
| Intercept | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Intercept[1] | 5.44*** | 1.57 | [ 3.06, 10.06] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 3521 | 2882 |
| Intercept[2] | 9.89*** | 3.14 | [ 5.45, 18.74] | 1.000 | [0.84, 1.20] | 0.000 | 1.002 | 3598 | 2951 |
| Intercept[3] | 23.10*** | 7.89 | [ 12.38, 46.79] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 3933 | 2899 |
| Intercept[4] | 81.88*** | 35.09 | [ 38.53, 199.63] | 1.000 | [0.84, 1.20] | 0.000 | 1.000 | 4569 | 3054 |
| Intercept[5] | 771.25*** | 625.85 | [202.52, 4855.18] | 1.000 | [0.84, 1.20] | 0.000 | 1.001 | 5246 | 2545 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing experienced | 1.31* | 0.14 | [ 1.04, 1.62] | 0.987 | [0.84, 1.20] | 0.214 | 1.001 | 3641 | 2884 |
| Daily pushing utilized (partner’s view) | 0.96 | 0.13 | [ 0.75, 1.25] | 0.615 | [0.84, 1.20] | 0.810 | 1.002 | 4044 | 2794 |
| Day | 1.60 | 0.73 | [ 0.68, 3.95] | 0.854 | [0.84, 1.20] | 0.182 | 1.002 | 4822 | 3245 |
| Own Actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Partner Actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Between-Person Effects | |||||||||
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing experienced | 2.32*** | 0.65 | [ 1.35, 3.99] | 1.000 | [0.84, 1.20] | 0.005 | 1.001 | 3394 | 2987 |
| Mean pushing utilized (partner’s view) | 1.25 | 0.35 | [ 0.70, 2.16] | 0.783 | [0.84, 1.20] | 0.365 | 1.001 | 3738 | 2562 |
| Random Effects | |||||||||
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Intercept) | 0.85 | 0.21 | [0.49, 1.34] | NA | NA | NA | 1.001 | 1569 | 2204 |
| sd(Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing experienced) | 0.17 | 0.14 | [0.01, 0.50] | NA | NA | NA | 1.004 | 1241 | 1651 |
| Additional Parameters | |||||||||
| sd(Daily pushing utilized (partner’s view)) | 0.16 | 0.14 | [0.01, 0.57] | NA | NA | NA | 1.003 | 1626 | 1714 |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| disc | 1.00 | 0.00 | [1.00, 1.00] | NA | NA | NA | NA | NA | NA |
plot(
bayestestR::p_direction(reactance_ordinal),
priors = TRUE
) +
coord_cartesian(xlim = c(-6, 6)) +
theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
## = 0.75), b_Intercept[4] and b_Intercept[3] (r = 0.82). This might lead
## to inappropriate results. See 'Details' in '?rope'.
Binary
introduce_binary_reactance <- function(data) {
data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
return(data)
}
df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
for (i in seq_along(implist)) {
implist[[i]] <- introduce_binary_reactance(implist[[i]])
}
}
formula <- bf(
is_reactance ~
pushing_self_cw + pushing_partner_cw +
pushing_self_cb + pushing_partner_cb +
#plan_self + plan_partner +
day +
# Random effects
(pushing_self_cw + pushing_partner_cw | coupleID)
#, autocor = autocor_str
)
prior1 <- c(
brms::set_prior("normal(0, 2.5)", class = "b")
, brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5)
, brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)
#brms::validate_prior(
# prior1,
# formula = formula,
# data = df_double,
# family = bernoulli()
# )
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]
is_reactance <- my_brm(
mi = use_mi,
imputed_data = implist,
formula = formula,
prior = prior1,
data = df_double,
family = brms::bernoulli(),
#control = list(adapt_delta = 0.95, max_treedepth = 15),
iter = iterations,
warmup = warmup,
chains = 4,
cores = 4,
seed = 7777,
file = file.path("models_cache_brms", paste0("is_reactance", suffix_onlypushing))
#, file_refit = 'always'
)## Warning: Rows containing NAs were excluded from the model.summarize_brms(
is_reactance,
stats_to_report = c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS'),
model_rows_fixed = model_rows_fixed,
model_rows_random = model_rows_random,
model_rownames_fixed = model_rownames_fixed,
model_rownames_random = model_rownames_random,
exponentiate = T) %>%
print_df(rows_to_pack = rows_to_pack)| OR | SE | 95% CI | pd | ROPE | inside ROPE | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|---|---|---|
| Intercept | 0.40*** | 0.11 | [0.23, 0.68] | 1.000 | [0.83, 1.20] | 0.002 | 1.000 | 3095 | 2830 |
| Within-Person Effects | |||||||||
| Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing experienced | 1.37* | 0.18 | [1.06, 1.87] | 0.993 | [0.83, 1.20] | 0.154 | 1.000 | 3238 | 2410 |
| Daily pushing utilized (partner’s view) | 1.16 | 0.26 | [0.80, 2.12] | 0.769 | [0.83, 1.20] | 0.509 | 1.000 | 2310 | 2170 |
| Day | 1.78 | 0.89 | [0.69, 4.81] | 0.874 | [0.83, 1.20] | 0.160 | 1.000 | 6984 | 2927 |
| Own Actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Partner Actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Between-Person Effects | |||||||||
| Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing experienced | 5.39*** | 2.58 | [2.38, 14.94] | 1.000 | [0.83, 1.20] | 0.000 | 1.001 | 2293 | 2628 |
| Mean pushing utilized (partner’s view) | 2.37 | 1.15 | [0.97, 6.82] | 0.972 | [0.83, 1.20] | 0.057 | 1.001 | 2556 | 2640 |
| Random Effects | |||||||||
| Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Intercept) | 1.52 | 0.31 | [1.03, 2.24] | NA | NA | NA | 1.001 | 1431 | 2494 |
| sd(Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing experienced) | 0.18 | 0.15 | [0.01, 0.56] | NA | NA | NA | 1.003 | 1256 | 1876 |
| Additional Parameters | |||||||||
| sd(Daily pushing utilized (partner’s view)) | 0.42 | 0.27 | [0.03, 1.15] | NA | NA | NA | 1.002 | 1292 | 1879 |
| sigma | NA | NA | NA | NA | NA | NA | NA | NA | NA |
plot(
bayestestR::p_direction(is_reactance),
priors = TRUE
) +
coord_cartesian(xlim = c(-6, 6)) +
theme_bw()## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length
conditional_spaghetti(
is_reactance,
effects = c('pushing_self_cw'),
group_var = 'coupleID',
plot_full_range = TRUE
)$pushing_self_cw
Report All Models
process_model_component <- function(obj) {
# Convert to string, modify, and evaluate
name <- deparse(substitute(obj))
if (report_hurdle) name <- paste0(name, '_hu')
if (report_ordinal) name <- paste0(name, '_ordinal')
return(get(name, envir = parent.frame()))
}
all_models <- report_side_by_side(
pa_sub,
pa_obj_log,
mood_gauss,
reactance_ordinal,
is_reactance,
stats_to_report = c('CI', 'pd'),
model_rows_random = process_model_component(model_rows_random),
model_rows_fixed = process_model_component(model_rows_fixed),
model_rownames_random = process_model_component(model_rownames_random),
model_rownames_fixed = process_model_component(model_rownames_fixed)
) [1] “pa_sub”
## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.[1] “pa_obj_log”
## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”
# pretty printing
summary_all_models <- all_models %>%
print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
add_header_above(
c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),
"Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5),
"Mood Gaussian" = (length(all_models) / 5),
"Reactance Ordinal" = (length(all_models) / 5),
"Reactance Dichotome" = (length(all_models) / 5))
)
export_xlsx(
summary_all_models,
rows_to_pack = process_model_component(rows_to_pack),
file.path("Output", paste0("AllModels", suffix_onlypushing, ".xlsx")),
merge_option = 'header',
simplify_2nd_row = TRUE,
line_above_rows = c(1,2),
line_below_rows = c(-1)
)
print(summary_all_models)| exp(Est.) pa_sub | 95% CI pa_sub | pd pa_sub | exp(Est.) pa_obj_log | 95% CI pa_obj_log | pd pa_obj_log | Est. mood_gauss | 95% CI mood_gauss | pd mood_gauss | OR reactance_ordinal | 95% CI reactance_ordinal | pd reactance_ordinal | OR is_reactance | 95% CI is_reactance | pd is_reactance | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Intercept | 57.28*** | [49.10, 66.97] | 1.000 | 124.24*** | [109.69, 139.58] | 1.000 | 3.78*** | [ 3.56, 4.00] | 1.000 | NA | NA | NA | 0.40*** | [0.23, 0.68] | 1.000 |
| Hurdle Intercept | 3.79*** | [ 2.27, 6.26] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Conditional Within-Person Effects | |||||||||||||||
| Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily pushing experienced | 0.97 | [ 0.91, 1.03] | 0.856 | 1.03 | [ 0.97, 1.08] | 0.807 | 0.02 | [-0.05, 0.08] | 0.670 | 1.31* | [ 1.04, 1.62] | 0.987 | 1.37* | [1.06, 1.87] | 0.993 |
| Daily pushing utilized (partner’s view) | 0.97 | [ 0.90, 1.03] | 0.847 | 1.02 | [ 0.97, 1.06] | 0.797 | 0.09* | [ 0.01, 0.17] | 0.984 | 0.96 | [ 0.75, 1.25] | 0.615 | 1.16 | [0.80, 2.12] | 0.769 |
| Day | 0.91 | [ 0.79, 1.06] | 0.894 | 0.97 | [ 0.87, 1.08] | 0.718 | 0.28*** | [ 0.11, 0.46] | 1.000 | 1.60 | [ 0.68, 3.95] | 0.854 | 1.78 | [0.69, 4.81] | 0.874 |
| Own actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Partner actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Daily weartime | NA | NA | NA | 1.00** | [ 1.00, 1.00] | 0.998 | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Conditional Between-Person Effects | |||||||||||||||
| Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Mean pushing experienced | 1.32* | [ 1.00, 1.71] | 0.977 | 0.97 | [ 0.86, 1.10] | 0.680 | 0.12 | [-0.08, 0.32] | 0.880 | 2.32*** | [ 1.35, 3.99] | 1.000 | 5.39*** | [2.38, 14.94] | 1.000 |
| Mean pushing utilized (partner’s view) | 1.18 | [ 0.90, 1.53] | 0.883 | 1.04 | [ 0.92, 1.18] | 0.739 | 0.08 | [-0.12, 0.30] | 0.786 | 1.25 | [ 0.70, 2.16] | 0.783 | 2.37 | [0.97, 6.82] | 0.972 |
| Mean weartime | NA | NA | NA | 1.00 | [ 1.00, 1.00] | 0.812 | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Within-Person Effects | |||||||||||||||
| Hu Daily persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pushing experienced | 1.55** | [ 1.17, 2.24] | 0.998 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily pushing utilized (partner’s view) | 1.84*** | [ 1.38, 2.68] | 1.000 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Day | 0.61 | [ 0.37, 1.00] | 0.974 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Own actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Partner actionplan | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Daily weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hurdle Between-Person Effects | |||||||||||||||
| Hu Mean persuasion experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean persuasion utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pressure experienced | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pressure utilized (partner’s view) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pushing experienced | 0.46** | [ 0.26, 0.81] | 0.997 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean pushing utilized (partner’s view) | 0.62 | [ 0.33, 1.08] | 0.952 | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Hu Mean weartime | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Random Effects | |||||||||||||||
| sd(Intercept) | 0.37 | [0.29, 0.50] | NA | 0.31 | [0.24, 0.42] | NA | 0.61 | [0.48, 0.79] | NA | 0.85 | [0.49, 1.34] | NA | 1.52 | [1.03, 2.24] | NA |
| sd(Hurdle Intercept) | 1.25 | [0.93, 1.74] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Daily pushing experienced) | 0.06 | [0.00, 0.16] | NA | 0.09 | [0.02, 0.17] | NA | 0.07 | [0.00, 0.16] | NA | 0.17 | [0.01, 0.50] | NA | 0.18 | [0.01, 0.56] | NA |
| sd(Daily pushing utilized (partner’s view)) | 0.09 | [0.01, 0.18] | NA | 0.03 | [0.00, 0.10] | NA | 0.14 | [0.06, 0.24] | NA | 0.16 | [0.01, 0.57] | NA | 0.42 | [0.03, 1.15] | NA |
| sd(Hu Daily persuasion experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily persuasion utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pressure experienced) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pressure utilized (partner’s view)) | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pushing experienced) | 0.46 | [0.07, 1.02] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| sd(Hu Daily pushing utilized (partner’s view)) | 0.51 | [0.10, 1.13] | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| Additional Parameters | |||||||||||||||
| sigma | 0.67 | [0.64, 0.71] | NA | 0.54 | [0.52, 0.57] | NA | 0.91 | [0.87, 0.94] | NA | NA | NA | NA | NA | NA | NA |
Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)
- beepr (version 2.0; Bååth R, 2024)
- R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
- R.oo (version 1.27.0; Bengtsson H, 2003)
- R.utils (version 2.12.3; Bengtsson H, 2023)
- brms (version 2.22.0; Bürkner P, 2017)
- digest (version 0.6.37; Eddelbuettel D, 2024)
- Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
- bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
- lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
- DHARMa (version 0.4.7; Hartig F, 2024)
- wbCorr (version 0.1.22; Küng P, 2023)
- see (version 0.9.0; Lüdecke D et al., 2021)
- tibble (version 3.2.1; Müller K, Wickham H, 2023)
- R (version 4.4.1; R Core Team, 2024)
- openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
- ggplot2 (version 3.5.1; Wickham H, 2016)
- forcats (version 1.0.0; Wickham H, 2023)
- stringr (version 1.5.1; Wickham H, 2023)
- rvest (version 1.0.4; Wickham H, 2024)
- tidyverse (version 2.0.0; Wickham H et al., 2019)
- readxl (version 1.4.3; Wickham H, Bryan J, 2023)
- dplyr (version 1.1.4; Wickham H et al., 2023)
- purrr (version 1.0.2; Wickham H, Henry L, 2023)
- readr (version 2.1.5; Wickham H et al., 2024)
- xml2 (version 1.3.6; Wickham H et al., 2023)
- tidyr (version 1.3.1; Wickham H et al., 2024)
- knitr (version 1.49; Xie Y, 2024)
- kableExtra (version 1.4.0; Zhu H, 2024)
- Bååth R (2024). beepr: Easily Play Notification Sounds on any Platform. R package version 2.0, https://CRAN.R-project.org/package=beepr.
- Bengtsson H (2003). “The R.oo package - Object-Oriented Programming with References Using Standard R Code.” In Hornik K, Leisch F, Zeileis A (eds.), Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003). https://www.r-project.org/conferences/DSC-2003/Proceedings/Bengtsson.pdf.
- Bengtsson H (2003). “The R.oo package - Object-Oriented Programming with References Using Standard R Code.” In Hornik K, Leisch F, Zeileis A (eds.), Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003). https://www.r-project.org/conferences/DSC-2003/Proceedings/Bengtsson.pdf.
- Bengtsson H (2023). R.utils: Various Programming Utilities. R package version 2.12.3, https://CRAN.R-project.org/package=R.utils.
- Bürkner P (2017). “brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software, 80(1), 1-28. doi:10.18637/jss.v080.i01 https://doi.org/10.18637/jss.v080.i01. Bürkner P (2018). “Advanced Bayesian Multilevel Modeling with the R Package brms.” The R Journal, 10(1), 395-411. doi:10.32614/RJ-2018-017 https://doi.org/10.32614/RJ-2018-017. Bürkner P (2021). “Bayesian Item Response Modeling in R with brms and Stan.” Journal of Statistical Software, 100(5), 1-54. doi:10.18637/jss.v100.i05 https://doi.org/10.18637/jss.v100.i05.
- Eddelbuettel D (2024). digest: Create Compact Hash Digests of R Objects. R package version 0.6.37, https://CRAN.R-project.org/package=digest.
- Eddelbuettel D, Francois R, Allaire J, Ushey K, Kou Q, Russell N, Ucar I, Bates D, Chambers J (2024). Rcpp: Seamless R and C++ Integration. R package version 1.0.13-1, https://dirk.eddelbuettel.com/code/rcpp.html, https://github.com/RcppCore/Rcpp, https://www.rcpp.org. Eddelbuettel D, François R (2011). “Rcpp: Seamless R and C++ Integration.” Journal of Statistical Software, 40(8), 1-18. doi:10.18637/jss.v040.i08 https://doi.org/10.18637/jss.v040.i08. Eddelbuettel D (2013). Seamless R and C++ Integration with Rcpp. Springer, New York. doi:10.1007/978-1-4614-6868-4 https://doi.org/10.1007/978-1-4614-6868-4, ISBN 978-1-4614-6867-7. Eddelbuettel D, Balamuta J (2018). “Extending R with C++: A Brief Introduction to Rcpp.” The American Statistician, 72(1), 28-36. doi:10.1080/00031305.2017.1375990 https://doi.org/10.1080/00031305.2017.1375990.
- Gabry J, Mahr T (2024). “bayesplot: Plotting for Bayesian Models.” R package version 1.11.1, https://mc-stan.org/bayesplot/. Gabry J, Simpson D, Vehtari A, Betancourt M, Gelman A (2019). “Visualization in Bayesian workflow.” J. R. Stat. Soc. A, 182, 389-402. doi:10.1111/rssa.12378 https://doi.org/10.1111/rssa.12378.
- Grolemund G, Wickham H (2011). “Dates and Times Made Easy with lubridate.” Journal of Statistical Software, 40(3), 1-25. https://www.jstatsoft.org/v40/i03/.
- Hartig F (2024). DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models. R package version 0.4.7, http://florianhartig.github.io/DHARMa/.
- Küng P (2023). wbCorr: Bivariate Within- and Between-Cluster Correlations. University of Zürich. R package version 0.1.22. https://github.com/Pascal-Kueng/wbCorr.
- Lüdecke D, Patil I, Ben-Shachar M, Wiernik B, Waggoner P, Makowski D (2021). “see: An R Package for Visualizing Statistical Models.” Journal of Open Source Software, 6(64), 3393. doi:10.21105/joss.03393 https://doi.org/10.21105/joss.03393.
- Müller K, Wickham H (2023). tibble: Simple Data Frames. R package version 3.2.1, https://CRAN.R-project.org/package=tibble.
- R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.
- Schauberger P, Walker A (2024). openxlsx: Read, Write and Edit xlsx Files. R package version 4.2.7.1, https://CRAN.R-project.org/package=openxlsx.
- Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.
- Wickham H (2023). forcats: Tools for Working with Categorical Variables (Factors). R package version 1.0.0, https://CRAN.R-project.org/package=forcats.
- Wickham H (2023). stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.5.1, https://CRAN.R-project.org/package=stringr.
- Wickham H (2024). rvest: Easily Harvest (Scrape) Web Pages. R package version 1.0.4, https://CRAN.R-project.org/package=rvest.
- Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). “Welcome to the tidyverse.” Journal of Open Source Software, 4(43), 1686. doi:10.21105/joss.01686 https://doi.org/10.21105/joss.01686.
- Wickham H, Bryan J (2023). readxl: Read Excel Files. R package version 1.4.3, https://CRAN.R-project.org/package=readxl.
- Wickham H, François R, Henry L, Müller K, Vaughan D (2023). dplyr: A Grammar of Data Manipulation. R package version 1.1.4, https://CRAN.R-project.org/package=dplyr.
- Wickham H, Henry L (2023). purrr: Functional Programming Tools. R package version 1.0.2, https://CRAN.R-project.org/package=purrr.
- Wickham H, Hester J, Bryan J (2024). readr: Read Rectangular Text Data. R package version 2.1.5, https://CRAN.R-project.org/package=readr.
- Wickham H, Hester J, Ooms J (2023). xml2: Parse XML. R package version 1.3.6, https://CRAN.R-project.org/package=xml2.
- Wickham H, Vaughan D, Girlich M (2024). tidyr: Tidy Messy Data. R package version 1.3.1, https://CRAN.R-project.org/package=tidyr.
- Xie Y (2024). knitr: A General-Purpose Package for Dynamic Report Generation in R. R package version 1.49, https://yihui.org/knitr/. Xie Y (2015). Dynamic Documents with R and knitr, 2nd edition. Chapman and Hall/CRC, Boca Raton, Florida. ISBN 978-1498716963, https://yihui.org/knitr/. Xie Y (2014). “knitr: A Comprehensive Tool for Reproducible Research in R.” In Stodden V, Leisch F, Peng RD (eds.), Implementing Reproducible Computational Research. Chapman and Hall/CRC. ISBN 978-1466561595.
- Zhu H (2024). kableExtra: Construct Complex Table with ‘kable’ and Pipe Syntax. R package version 1.4.0, https://CRAN.R-project.org/package=kableExtra.